## Examine spud to production time distributions ----
dt[, .N, by=.(is.na(spud_to_prod), is.na(spud_date))]
dt[, .N, by=sign(spud_to_prod)]

StP.cutoff.on = 24 # 2 year primary lease terms
StP.cutoff.off = 10*12 # 10 year leases

dt[offshore==0, reasonable := 1*between(spud_to_prod, 0, StP.cutoff.on)]
dt[offshore==1, reasonable := 1*between(spud_to_prod, 0, StP.cutoff.off)]

# Basically, the reasonable ones have higher production, holding constant fed.land-x-offshore.
# The missing ones aren't really producing anymore?
palette(alpha(col_lookup, 1)) # Order: gas fed, gas nonfed, oil fed, oil nonfed
plot(1:4, pch=19, cex=5, col=1:4)

dt.plot = dt[reasonable==1, .(fed.land, offshore, prod_type, spud_to_prod)]

dev.off()
bw=3
i=0
for (i in 0:1) { # onshore/offshore
  pdf(output.dir%&%'/Spud_to_Prod_Dists_'%&%ifelse(i==0, 'Onshore','Offshore')%&%'_'%&%today()%&%'.pdf', family='serif', width=11, height=8.5)
  dt.plot[offshore==i & fed.land==1 & prod_type=='Gas'][,  plot(density(spud_to_prod, bw=bw, from=0, to=ifelse(i==0, StP.cutoff.on, StP.cutoff.off)), 
                                                                main='', xlab='Months from Spud to Production', cex=1.5, col=1, lwd=3, ylim=c(0, 0.12),
                                                                cex.lab=1.5, cex.axis=1.5, xaxt=ifelse(i==0, 'n', 's'))]
  grid(nx=NA, ny=NULL) # draw y grid
  if (i==0) { axis(1, at=seq(0, StP.cutoff.on, by=6), cex.axis=1.5); abline(v=seq(0, StP.cutoff.on, by=6), lty=3, col='lightgray')
  } else grid(nx=NULL, ny=NA) # if we didn't do it manually, draw x grid
  dt.plot[offshore==i & fed.land==0 & prod_type=='Gas'][, lines(density(spud_to_prod, bw=bw, from=0, to=ifelse(i==0, StP.cutoff.on, StP.cutoff.off)), col=2, lwd=3)]
  dt.plot[offshore==i & fed.land==1 & prod_type=='Oil'][, lines(density(spud_to_prod, bw=bw, from=0, to=ifelse(i==0, StP.cutoff.on, StP.cutoff.off)), col=3, lwd=3)]
  dt.plot[offshore==i & fed.land==0 & prod_type=='Oil'][, lines(density(spud_to_prod, bw=bw, from=0, to=ifelse(i==0, StP.cutoff.on, StP.cutoff.off)), col=4, lwd=3)]
  
  legend('top', legend=c('Oil Wells, Nonfederal',
                         'Oil Wells, Federal',
                         'Gas Wells, Nonfederal',
                         'Gas Wells, Federal'),
         col=c(4,3,2,1), lty=1, lwd=3, cex=1.5, bg='white')
  dev.off()
}
rm(dt.plot)
rm(states_sp)

# Save histograms/densities
stp.densities <- stp.densities.2018 <-  array(NA, dim=c(StP.cutoff.off, 2, 2, 2),
                                              dimnames=list(1:dt[reasonable==1, max(spud_to_prod)],
                                                            c('Oil','Gas'),
                                                            c('Nonfederal','Federal'),
                                                            c('Onshore','Offshore')))
for (p in 1:2) {
  for (f in 1:2) {
    for (o in 1:2) {
      h = dt[reasonable==1][prod_type==dimnames(stp.densities)[[2]][p] &
                              fed.land==f-1 & offshore==o-1, 
                            hist(spud_to_prod, breaks=0:StP.cutoff.off)]
      h.2018 = dt[first_prod_year==2018][reasonable==1][prod_type==dimnames(stp.densities)[[2]][p] &
                                                          fed.land==f-1 & offshore==o-1, 
                                                        hist(spud_to_prod, breaks=0:StP.cutoff.off)]
      stp.densities[,p,f,o] = h$density
      stp.densities.2018[,p,f,o] = h.2018$density
    }
    
  }
}
dimnames(stp.densities.2018)
# Check that they sum to 1
apply(stp.densities, 2:4, sum) 
apply(stp.densities.2018, 2:4, sum)
